Genome-wide identification and characterization of the HD-Zip gene family and expression analysis in response to stress in Rehmannia glutinosa Libosch

ABSTRACT The HD-Zip family of transcription factors is unique to the plant kingdom, and play roles in modulation of plant growth and response to environmental stresses. R. glutinosa is an important Chinese medicinal material. Its yield and quality are susceptible to various stresses. The HD-Zip transcription factors is unique to the plant, and roles in modulation of plant growth and response to environmental stresses. However, there is no relevant research on the HD-ZIP of R. glutinosa. In this study, 92 HD-Zip transcription factors were identified in R. glutinosa, and denominated as RgHDZ1-RgHDZ92. Members of RgHDZ were classified into four groups (HD-ZipI-IV) based on the phylogenetic relationship of Arabidopsis HD-Zip proteins, and each group contains 38, 18, 17, and 19 members, respectively. Expression analyses of RgHDZ genes based on transcriptome data showed that the expression of these genes could be induced by the endophytic fungus of R. glutinosa. Additionally, we showed that RgHDZ genes were differentially expressed in response to drought, waterlogging, temperature, and salinity treatments. This study provides important information for different expression patterns of stress-responsive HD-Zip and may contribute to the better understanding of the different responses of plants to biotic and abiotic stresses, and provide a molecular basis for the cultivation of resistant varieties of R. glutinosa.


Introduction
R. glutinosa and its active principles possess wide pharmacological actions on the blood system, immune system, endocrine system, cardiovascular system, and the nervous system. R. glutinosa yield and quality are limited by biotic, such as root rot, wheel spot, and fusarium wilt caused by microorganisms, as well as abiotic stresses, such as salt, heavy drought, and waterlogging. Exploring the molecular foundation of the generegulatory systems underlying plant responses to both abiotic and biotic stresses is crucial for medicinal plants improvement. In this regard, homeodomain leucine zipper (HD-Zip) genes family concerned with enlightening plant growth and tolerance to environmental stresses are considered key players for plants improvement. 1 Therefore, excavation stress-related genes can provide a theoretical basis for breeding R. glutinosa resistant varieties. Ma Ligang et al. reported a reference genome of R. glutinosa using Nanopore technology, Illumina, and Hi-C sequencing. The assembly genome is 2.49 Gb long with a scaffold N50 length of 70 Mb and high heterozygosity (2%). 2 This laid the foundation for the discovery of genes related to R. glutinosa stress. Transcription factors can specifically bind to specific sequences of target genes, activate, or inhibit the expression of target genes, and then regulate a variety of life processes. The HD-Zip transcription factors is unique to the plant, and roles in modulation of plant growth and response to environmental stresses. 3 It has been identified in a variety of plants, such as 33 in grapes (Vitis vinifera L.), 4 55 in corn, 5 48 in Arabidopsis thaliana, 5 45 in sesame, 3

and 63 in
Populus trichocarpa. 6 HD-ZIP family transcription factors have two domains, homeodomian (HD) and Leucine-Zip (LZ). The HD is 60 amino acids in length and adopts a structure of three αhelices connected by a loop and a turn. The leucine-zipper is formed by five or more highly conserved leucine residues and, in general, exactly located every six amino acids. The HD-Zip family further divided into four subfamilies, HD-ZipI-IV, according to the conserved HD-Zip domain, additional conserved domain, structure features, and functions. HD-Zip class I protein has the simplest structure, which is closely related to plant abiotic stress, ABA response, blue light, de-etiolated, and embryogenesis. 7,8 Overexpression of MdHB-7 promoted endogenous ABA accumulation, stomatal closure, and reactive oxygen species (ROS) detoxification in response to drought, whereas suppression of MdHB-7 expression had the opposite effect. 9 A quantitative real-time PCR analysis indicated Phehdz1, a HD-Zip of Phyllostachys edulis expression was significantly induced by drought, high salinity, and abscisic acid (ABA). 10 HD-ZipII TFs have been documented to affect plant growth and development mainly by regulating auxin, and have a distinguishing feature in their C-terminus, the CPSCE motif responsible for redox regulation of protein activity. [11][12][13][14] In addition, the HD-ZIPII genes also responds to biotic and abiotic stress in plants. 15,16 HD-ZipIII proteins can be distinguished from HD-ZipIV proteins by the presence of C-terminal MEKHLA domain which is absent in the HD-ZipIV proteins. 17 HD-ZipIII TFs are mainly involved in the formation of plant vascular bundles, auxin transport, 18 and organ development. [19][20][21] In addition, they are also involved in plant abiotic stress tolerance, such as OsHB4 modulates cadmium tolerance and accumulation in rice. HD-ZipⅣ TFs play vital role in the development of plant epidermis. 22, 23 Cai found that the HD-ZipIV transcription factor GL2-LIKE regulates male flowering time and fertility in cucumber. 24 In addition to regulating epidermal development, HD-ZipIV is also involved in plant abiotic stress response. 25 In recent years, miRNAs have been noticeable with their important biological functions. MicroRNAs (miRNAs), a class of endogenous small regulatory non-coding RNAs, are usually 20-24 nucleotides long that play pivotal roles in gene regulation during plant development. 26,27 Studies have shown that HD-ZipIII subfamily TFs are regulated by miRNAs to affect plant growth and development and resisitance. 28,29,30 HD-ZipIII genes are the main target genes of miR166, which have been down-regulated by miR166 to regulate the response of tomato seedlings to cold stress. 31 MiR166 plays a role in biotic and abiotic stress responses by regulating the transcription level of HD-ZipIII in Arabidopsis. 32 miRNAs are able to cleave their target mRNA of HD-Zip genes, thus regulating the functions of these genes. Therefore, predicting the miRNA that regulates HD-Zip genes is the important method for studying biological functions of HD-Zip TFs.
At present, there have been no relevant research on the HD-Zip transcription factors of R. glutinosa. In this study, R. glutinosa was used as the material to identify and explore the functions of HD-Zip genes. It provided a theoretical basis for revealing the anti-adversity regulation mechanism of R. glutinosa, which through qRT-PCR analysis, we explored the response of RgHDZ to different abiotic stresses, deepened our understanding of RgHDZ.

Physicochemical properties and structure analysis of R. glutinosa HD-Zip transcription factor
According to HMMER search, a total of 92 non-redundant R. glutinosa HD-Zip encoding genes were identified denominated as RgHDZ1-RgHDZ92. RgHDZs encode proteins of 163-858 amino acids, and the molecular weight of RgHDZ is 18.99kDa-93.91 kDa. Except for Class II HD-Zips TFs, the isoelectric points of most proteins of other subfamily are less than 7. It was indicated that Class II is composed of more basic amino acids, while other proteins contain more acidic amino acids. In addition, the average hydrophobicity of RgHDZ proteins is less than 0, indicating that they are all hydrophilic proteins (Table 1). Subcellular localization analysis revealed that RgHDZ proteins are located in the nucleus (Supplementary Table 1). Transmembrane analysis showed that RgHDZ proteins are non-transmembrane protein.

Phylogenetic analysis amongst the R. glutinosa
To reveal the phylogenetic relationships among the RgHDZ proteins, an unrooted phylogenetic tree was created to assess the genetic relationships between Arabidopsis and R. glutinosa HD-Zip proteins. The sequences of HD-Zip proteins of Arabidopsis were obtained from NCBI. As shown in Figure 1, these proteins can be divided into four distinct groups (HD-ZipI-IV), which is similar to that described in previous studies. Number of HD-ZipI-IV members in R. glutinosawas 38, 18, 17 and 19 respectively ( Figure 1). The results provide an important basis for functional prediction of HD-Zip proteins in R. glutinosa.

Motifs and chromosome localization analysis of RgHDZs
Using RgHDZ proteins sequences to construct a phylogenetic tree to study the evolutionary relationship among members of the family, the results are consistent with the phylogenetic analysis. According to the results of annotation, the structure map of RgHDZ genes was constructed. The results showed that the number of exons of RgHDZ was among 2-19 and the number of introns was between 1 and 18. And we found that subfamily III has the largest number of exons and introns, the number of exons is [16][17][18], and the number of introns is 16-19 ( Figure 2). To further study the origin and evolutionary pattern of RgHDZs, their amino acid sequences were subjected to the MEME tool and a total of 12 conserved motifs were identified. The identified motifs ranged from 6 to 200 amino acids in length motif 1 and 2, corresponding to the homeobox domain, and motif 5, corresponding to the LZ domain, were in common among all of the RgHDZs (SupplementaryFigure 2). Motifs 3, 4 and 6, which correspond to START domain, were shared in group III and group IV, while motif 8, which correspond to MEKHLA domain, were present in group III but absent in group IV.
(A) Glutinosa is an autotetraploid (4 n = 56), the difference between each set of chromosomes is very small and it is difficult to distinguish the two sets of chromosomes using Hi-C. Hence, only one set of the genome size was mounted to the chromosome level. 2 The location information shows that the RgHDZ gene family is unevenly distributed on the chromosomes of R. glutinosa. There are 43 genes distributed on 13 chromosomes. The number of RgHDZ genes on chromosomes 6 and 7 is the largest, with 14 in total, while only one RgHDZ gene distributed on chromosomes 8, 11, 12 and 13.

Cis-acting element prediction of RgHDZ
The promoter sequence of the RgHDZ genes (1.5 kb upstream of 5ʹUTR) was submitted to the PlantCARE database for ciselements prediction, and 72 cis-acting elements were identified. As shown in Supplementary figure 2, all RgHDZ genes contain 2 conventional promoter elements (TATA-box and CAAT-box). The cis-acting elements frequently predicted in most RgHDZ genes mainly fall into four categories: light-responsive elements, shedding Acid response element, anaerobic inducing element and MeJA response element. Among them, there are up to 24 types of light-responsive elements (3-AF1 binding site, AAAC-motif, ACA-motif, ACE, etc.). In addition, there were some that appear less frequently and were only predicted in a few RgHDZ genes, such as gibberellin response (GARE-motif, P-box, and TATA-box), salicylic acid response (TCA, SARE), auxin Elements (AuxRE, AuxRR-core) , circadian elements (circadian), endosperm expression (AACA-motif, GCN4-motif), and a cis-element (MBSI) is involved in the MYB binding of flavonoid biosynthesis gene regulation Site.

Synteny analysis of RgHDZ genes
In order to clarify the underlying evolutionary mechanism of the R. glutinosa RgHDZ gene family, we studied its replication events and the results showed that some large fragments of repetitions were detected on the chromosomes of R. glutinosa. These homologous gene pairs were   In order to clarify the role of selection pressure in the evolutionary species of RgHDZ family, the synonymous substitutions (Ks), the nonsynonymous substitution rate (Ka) and their ratio (Ka/Ks) of RgHDZ homologous genes were analyzed using TBtools. Results The Ka/Ks of the 19 pairs of homologous genes were all far less than 1 ( Table 2), indicating that the RgHDZ family has undergone strong purification selection during the evolution process

Prediction of protein-protein interaction network
In order to further understand the interaction of pepper RgHDZ proteins, an interactive network based on Arabidopsis orthodoxy was established using STRING. As shown in Figure 3 (ENHANCER OF TRY AND CPC 1). In interacting proteins networks, the largest hub is KAN proteins (Figure 3). Surprisingly, all class III HD-Zip of R. glutinosa interact with KAN. The class III HD-Zip and KANADI genes comprise a genetic system that patterns abaxial-adaxial polarity in lateral organs produced from the apical meristem. 33 Besides, we found that GL3, TTG1, EGL3, MYB66, CPC, TRY, MYB23, and ETC1 are all related to the development of plant epidermis. Such as GL3 and EGL3 participate in an intercellular regulatory circuit that controls cell patterning in the Arabidopsis root epidermis. 34 In addition, MYB66 and CPC not only participate in the development of Arabidopsis trichomes, but also regulate the biosynthesis of anthocyanins. [35][36][37][38] We found that there is a close interaction between II and III class HD-Zip proteins, which may be related to the fact that HD-ZipII and III binding sites share the same core sequence [AAT(G/C)ATT], thus suggesting that members of the two protein families may regulate common target genes. 39 It is speculated that RgHDZ is involved in the growth and development of R. glutinosa and the process of adversity stress through interaction protein analysis.

Analysis of miRNA regulating RgHDZ genes
To further analyze the biological functions of RgHDZ,we predicted the miRNA that may regulate the HD-Zip genes of R. glutinosa. Results a total of 498 miRNA were predicted to be distributed in 93 families. We found that the class I HD-Zip genes are mainly regulated by miR169 and miR159. Research has shown that the involvement of ma-miR169a (Musa acuminata L.) and ma-miR169b in the banana response to Foc4 (Fusarium oxysporumf. sp. cubense tropical race 4). 40 Under drought stress, the differential expression of miRNA regulates the expression of their target genes, resulting in multiple responses of physiological and biochemical pathways relative to drought tolerance of maize. 41 Several studies have shown that miR169, miR159 had the contributions to plant response toward drought, salinity, cold, and heat, respectively. 42 This corresponds to the biological function of the aforementioned HD-ZipI subfamily in response to abiotic stress. After prediction, it was found that there are 34 miRNA that regulate HD-ZipII genes, mainly distributed in 10 families: miR4221, miR6459, miR156, miR1072, miR169, miR5139, miR159. The miR6459 is a female characteristic expression gene. 43 The alterations of miR5139, contributed to the attenuation of Cd translocation into the shoot of QLQ. 44 In addition, we found that miR169 and miR159 family genes related to plant abiotic stress. The most number of miRNA were predicted of class III HD-Zip, 272 are distributed in 45 families, including miR166, miR165, miR5168, miR8029, miR447, miR172, etc. The most important ones are miR166 and miR165 families with similar biological functions. Studies have shown that the main target genes of miR166 and miR165 is HD-ZipIII, which regulates the response of target genes to plant stress through differential expression. 41,45 It has been predicted  that 101 miRNA regulating HD-ZipIV genes are distributed in 34 families, mainly miR9671, miR5651, miR7731, miR1149.2, miR9735, miR6169, miR7520, miR2619, miR5260, miR6222, miR482. Except for cellulose hydrolysis, all other cell wall degradation modules are predicted to be targeted by miR6222. 46 SlymiR482e-3p mediates tomato wilt disease by modulating ethylene response pathway. 47 It is known that miRNA can inhibit translation or degradation through binding sites (miRNA response elements; MREs) and mRNA complementary pairing. 48 As shown in Figure 4 -b, miR166-3p and HD-ZipIII gene sequences is reversed at this position complementary, indicating that R. glutinosaHD-ZipIII genes can be translationally inhibited or degraded by miR166-3p.
Through the combined analysis of R. glutinosa miRNA (PRJNA781675), transcriptome (PRJNA781424), and degradation (PRJNA785421) sequencing. We found that RgHDZ69 was cleaved by miR166a (Figure 4 -c) during the infection process of R. glutinosa endophytic fungus GG22 to participate in the regulation of the biological stress process ofR. glutinosa. Moreover, we found that the change trend of miR166 expression and RgHDZ69 were mutually antagonistic after GG22 infection (Figure 4 -d; Figure 5).

Expression pattern analysis of RgHDZ in response to endophytic fungi GG22 infection
In order to explore the response of R. glutinosa HD-Zip gene in biotic stress. Tissue culture seedlings of R. glutinosa that was co-cultured with endophytic fungi GG22 for 3, 12, 29 days, another group was not cocultured with GG22 were used as materials, named GG3, GG12, GG29 and CK3, CK12, CK29, respectively, for transcriptome sequencing. Set three biological replicates for each group, and analyze and process the data obtained after sequencing.
According to the transcriptome datas, the expression levels of the 92 RgHDZ genes were, generally, up or down-regulated expression after infection by the endophytic fungi GG22 ( Figure 5). According to statistics, the expression of the 50/92 RgHDZ genes were down-regulated during the three sampling periods, and the performance of the class III and IV RgHDZ genes were particularly obvious. The differential expression of the group indicates that the HD-Zip family genes may be involved in the regulation of the stress response of R. glutinosa. The functions of homologous genes are known to be similar. In order to explore more different functions of RgHDZ genes, non-homologous genes were selected for experiments according to the results of phylogenetic tree and homologous gene analysis. For the purpose of studying whether it is involved in the abiotic stress response of R. glutinosa, 11 differentially expressed genes were selected from different branches of each subfamily (Figure 1), and qRT-PCR was used to detect their effects under drought, waterlogging, high temperature, low temperature, and salt stress to analyze the expression pattern of RgHDZ genes in response to abiotic stresses. Figure 6. Expression of RgHDZgenes in response to high and low temperature stresses. Total RNAs from the stressed samples were extracted from the third leaf of each plant. All data represent the mean ± standard deviation of three independent experiments. abcde indicate significant differences than control group (Student's t-test; P < .05).

Expression pattern analysis of RgHDZgenes in response to high and low temperature stresses
The expression levels of the selected RgHDZ genes were changed in response to high and low-temperature stresses, but some differences were present among these genes. Under hightemperature stress, the expression levels of RgHDZ45 and RgHDZ65 were significantly up-regulated through all the time points, while RgHDZ45 and 65 were up-regulated through all the time points in low-temperature stresses, with the highest expression levels at 48 h. In addition, some RgHDZ genes exhibited similar expression profiles under high-and lowtemperature stresses, such as RgHDZ58, RgHDZ82and RgHDZ92. However, some RgHDZ genes displayed an obvious decrease in expression under high-and low-temperature stresses at certain time points, such as RgHDZ7 and RgHDZ20 at 12 h and 24 h, respectively ( Figure 6). Expression of RgHDZ genes in response to salinity stresses. Total RNAs from the stressed samples were extracted from the third leaf of each plant. All data represent the mean ± standard deviation of three independent experiments. abcde indicate significant differences than control group (Student's t-test; P < .05). e2096787-10 Y. ZHU ET AL.

Expression pattern analysis of RgHDZ genes in response to salinity stresses
In order to identify the role of RgHDZ genes under salt stress, the expression patterns of the above 11 genes that respond to temperature stress were analyzed their expression patterns by qRT-PCR. Overall, the expression levels of all the selected genes were significantly changed in response to salinity stress, but some differences were present among these genes. Under salinity stress, the expression levels of  All data represent the mean ± standard deviation of three independent experiments. abcde indicate significant differences than control group (Student's t-test; P < .05).

Expression pattern analysis of RgHDZ genes in response to drought and waterlogging stresses
It was observed that most of the RgHDZ genes were differentially affected under drought and waterlogging stresses. For example, two RgHDZ genes, RgHDZ7 and RgHDZ45, were upregulated under drought stress, with their expression markedly increasing at 12 h, and peaking at 48 h (Figure11). RgHDZ20 was up-regulated under waterlogging stress, with their expression markedly increasing at 12 h to 48 h, and peaking at 12 h. Some RgHDZ genes, such as RgHDZ20, RgHDZ58, RgHDZ88, and RgHDZ92, were also up-regulated under stress, but their expression levels were peaked at 12 h, implying their roles in early response to drought and Waterlogging stress. However, some RgHDZ genes, such as RgHDZ20, RgHDZ26, RgHDZ45, Figure 9. Expression of RgHDZ genes in response to Fusarium solani stresses. Total RNAs from the stressed samples were extracted from the third leaf of each plant. All data represent the mean ± standard deviation of three independent experiments. abcde indicate significant differences than control group (Student's t-test; P < .05).
RgHDZ65, RgHDZ58, RgHDZ88, RgHDZ82, and RgHDZ92, exhibited similar expression profiles after drought and Waterlogging stresses. Moreover, it was showed that the response of RgHDZ genes to waterlogging is lagging to drought, such as RgHDZ7 and RgHDZ67 (Figure 8).

Expression pattern analysis of RgHDZ genes in response to biotic stresses
Different from GG22, Fusarium solani, an important plant pathogen can cause root rot in plants, which mainly infects the rhizomes of the host, causing root rot and blocking the ducts, and ultimately the host death. In this experiment, the change of RgHDZgenes expression in the tissue cultured seedlings of R. glutinosa co-cultured with Fusarium solani was investigated to study the role of RgHDZ genes in fungi stress of R. glutinosa. Expression of RgHDZ genes change significantly under pathogenic fungus stress, but it is generally at a low level (Figure 9). Overall, the expression levels of RgHDZ genes were significantly down-regulated in response to pathogenic fungus stress, but RgHDZ20, RgHDZ45, RgHDZ82, and RgHDZ88 up-regulated at certain time points, with the highest expression levels at 3 d, 3 d, 12 d, and 12 d, respectively. It can also be found from the results that the class III HD-Zip genes RgHDZ65, RgHDZ58, RgHDZ67, and RgHDZ71 were down-regulated in all sampling periods after Fusarium solani infection, which is similar with nonpathogenic fungus GG22 infection. Moreover, RgHDZ7 and RgHDZ20 exhibited similar expression profiles after pathogenic fungus stress. Taken together, these results suggested that these genes might play a vital role in response to multiple biotic stresses in R. glutinosa (Figure 9).
↑ means at least one time of expression that is up-regulated and there is no down-regulation in any other periods, ↓means at least one time of expression that is down-regulated and there is no up-regulation in any other periods, M means both upand down-regulation of expression in different periods,means no significant change in expression in all periods.
As seen in Table 3, all 11 RgHDZ genes were able to respond to a variety of different stress conditions. In general, the majority of RgHDZ gene expressions were up-regulated in response to abiotic stresses like heat, salt, and drought; however, the majority of RgHDZ gene expressions exhibited a decreasing trend in response to pathogenic fungal infestation.

Identification and phylogenetic analysis of HD-Zip TFs in R. glutinosa
HD-Zip genes encode a family of plant-specific transcription factors involved in various biological processes in plants. In this study, we identified 92 HD-Zip genes from R. glutinosa at the Genomic level. The number of HD-Zip members in R. glutinosa was nearly twice as much to that of Arabidopsis (48), 17 rice(48), 49 wheat (46). 25 Combined with the analysis of collinearity, it was found that it may be due to the replication of large fragments that caused the amplification of the HD-Zip genes family in R. glutinosa.
At present, more and more evidences show that the HD-Zip genes family plays an important role in plant life. For example, AtHB13 was able to confer tolerance to freezing temperatures via the induction of glucanase (GLU and PR2) and chitinase (PR4) proteins, 5 which could also regulate drought stress together with the target JUB1 gene. 50,51 ATHB23 plays role in salt and osmotic stress tolerance in Arabidopsis. 52 ATHB16 affects the content of ABA in plants by responding to blue light, 52 ATHB4 and ATHAT3, members of the class II HD-Zip transcription factor, have been shown to play an instrumental role in the responses to shade. 53 However, REV in Arabidopsis is related to the development of leaf polarity and overlaps with PHB and PHV functions. 54 Many studies have shown that class IV HD-Zip proteins are involved in plant abiotic stress processes, such as overexpression of AtHDG11 gene confers drought tolerance in plant; 55,56 AtPDF2, the Arabidopsis defensin gene, promotes cytoplasmic Cd efflux via chelation, thereby enhancing Cd detoxification and apoplastic accumulation; 57 According to the phylogenetic tree, collinearity analysis and Ka/Ks analysis, we can speculate that the RgHDZ gene is involved in regulating the growth and development of R. glutinosa and adversity stress.

RgHDZs protein interaction network
The establishment of protein interaction network is an important method for predicting protein function, which can be predicted based on the protein directly or indirectly interacting with the protein to be tested. As shown in Figure 5, there were interactions between RgHDZs protein members, which jointly regulate the complex life activities of plants. In addition, many proteins that do not belong to the HD-Zip interact with RgHDZs, such as KAN, GL3, TTG1, EGL3, MYB66, CPC, TRY, MYB23, and ETC1. The core protein KAN in the interaction network diagram is related to the development of plant organs, 33,58 and most of the proteins that interact with KAN belong to the class III RgHDZ protein. According to the literature search, we know that the class III HD-Zip protein is related to the development of plant epidermis, so through interaction analysis, not only helps us understand the protein interaction relationship of RgHDZ, but also can further speculate on its biological function. For example, we can speculate that type III RgHDZ may play an important role in the development of the organs of R. glutinosa. Based on the proven biological functions of these proteins, we predict the biological function of the RgHDZs protein to be tested, which may respond to a variety of abiotic stress processes, such as controlling the ABA content in R. glutinosa to affect plant resistance ability.

MicroRNA prediction
MiRNA, short 21-nucleotide RNA molecules, play an important role in post-transcriptional regulation of gene expression miRNA prediction of RgHDZ genes found that different HD-Zip family genes will be regulated by the same miRNA, and the same target gene will also be regulated by different miRNA. At present, there were more studies on the targeted regulation of miR166/miR165 on HD-ZipIII family genes. The Arabidopsis miR166/165 group targets five members of the HD-ZipIII transcription factor (REVOLUTA, PHABULOSA, PHAVOLUTA, CORONA, ATHB15, and ATHB8) functions by cleaving target mRNAs through complementary base pairing. 59,60,61 Furthermore, studies have shown that the differential expression of miR166/165 could regulate the abiotic response of plants. 45 Kitazumi researched that Inverse relationships have been established between the expression of two salt stress-regulated miRNA (miR166, miR159) and the transcriptional regulators they control (HD-Zip-Phabulosa/Phavoluta and Myb101, respectively) . Moreover, the expression of miR166, miR159 and miR156 was upregulated in Cavendish banana roots inoculated with Fusarium oxysporum. 62 Apparently, miR166 participates in the regulation of a variety of life activities, and also plays an important role in the process of plant stress response. According to the prediction results, we have found that the HD-ZipIII subfamily of R. glutinosa was regulated by miR166 and miR165, indicating that HD-ZipIIIs participates in the biotic and abiotic stress process of R. glutinosa and is regulated by miR166 and miR165. In addition, the other three subfamily HD-Zip genes have predicted miRNA related to plant stress. All in all, the RgHDZ genes may be regulated by miRNA in the process of participating in the stress response.

The expression pattern of RgHDZ genes
The expression of RgHDZ genes under biotic and abiotic stresses were tested, and the results clearly verified the hypothesis that it participates in the stress of R. glutinosa. (Table 3) Observing the expression of 11 RgHDZ genes under different abiotic stresses, it was found that most of the genes after abiotic stress were up-regulated compared to the control group during most of the period. And Compared with temperature and salt stresses, the RgHDZ genes had the most significant up-regulation reaction occurred in the drought and waterlogging stress group. This indicates that RgHDZ genes may play an important role in the process of R. glutinosa coping with drought and waterlogging stress. Related studies in Arabidopsis also found that the expression of HD-ZipI subfamily genes ATHB6 and ATHB7 were induced by water stress, 63,64 there were few studies on water-related stress in the III and IV subfamily. Our studies have found that its expression is induced in water-related stress, indicating that it may be involved in abiotic stress related to water in R. glutinosa.
Different from abiotic stress expression pattern, that the changes of RgHDZ genes expression after infection by Fusarium solani was very similar to infection by GG22, compared with the control group, the expression of RgHDZ was mainly down-regulated. Among them, RgHDZ88 and RgHDZ82 were down-regulated in the three sampling times. Therefore, the RgHDZ genes of R. glutinosa may mainly respond to the biological stress of fungal infection by down-regulating its expression.

Phylogenetic, exon-intron structure and conserved motif analyses of HD-zip family in R. glutinosa
The protein sequences of HD-Zip from R. glutinosa and Arabidopsis were used to construct the phylogenetic tree by MEGA 7, using the neighbor-joining (NJ) method with 1000 bootstrap replications. Import the results into online software iTOL (https://itol.embl.de/) for beautification. Use GSDS 2.0 (Gene Structure Display Server) (http://gsds.gaolab.org/) to determine the intron-exon structure of the RgHDZ gene. Conserved motifs present in RgHDZ were identified using MEME (https://meme-suite.org/meme/tools/ meme). The following parameters were used: number of repetitions any, maximum number of motifs 12, and the optimum motif widths were constrained to between 6 and 100 residues.

Promoter analysis
The 1500 bp upstream sequences of coding region of RgHDZ genes were downloaded from R. glutinosa Genome Database. The cis-regulatory elements were identified using online program PlantCARE(http://bioinformatics.psb.ugent.be/webt ools/plantcare/html/).

Interaction protein analysis
Interaction protein in RgHDZs were identified using STRING (https://string-db.org/), and used Cytoscape to draw a protein interaction network diagram.
Six groups of R. glutinosa tissue cultured seedlings RNA were extracted, and after passing the quality test, the sRNA library was constructed, using the Hiseq 2500 platform for sequencing. Then we processed and analyzed the raw data of sequencing, removed the 3ʹlinker and junk sequence, and obtained clean data; The sequences with a base length of 18-25 nt are reserved for comparison with various RNA databases, such as mRNA, RFam and Repbase database, and filtered. The finally obtained data is valid data for subsequent small RNA data analysis. Compare the valid data with the miRNAs included in miRBase for miRNAs identification.
Perform target gene prediction on the identified miRNAs, and verify the predicted target genes by degradome sequencing. Construction of degradation group library: (1) Capture mRNA with magnetic beads and connect with 3' and 5' adaptors, (2) Mixed reverse transcription of Biotinylated Random Primers and mRNA, (3) PCR amplification, after completing the entire library preparation, the constructed library is sequenced with Illumina Hiseq2000/2500, and the sequencing read length is single-ended 1 × 50 bp.
The target gene corresponding to the predicted miRNA is combined with the mRNA in the generated degradation group density file to find the common mRNA, which is the target gene of the miRNA. The peak classification and score of the degradation group are given, and the predicted results are plotted (t-plots).

Experimental materials
Fusarium solani:In the early stage, the laboratory separated and identified the fungi species from the leaves of R. glutinosa, PDA culture medium 5 to 7 days activation standby.
Transfer the R. glutinosa tissue culture seedlings grown in the MS medium to a plant height of 3-5 cm to the vermiculite medium to adapt to grow for 3 days, and then subjected to high (40°C), low (0°C) temperature, drought, waterlogging, and salt stress (100 mM, 200 mM). In addition, the biological stress group was inoculated with a 5 mm × 5 mm piece of Fusarium solani. Three biological replicates were performed for each experiment, and the sample of each replicate from 15 shares plants.

Expression characteristics of RgHDZ after GG22 infection
To analyze the expression profiles of RgHDZ genes in response to GG22 (R. glutinosa endophytic fungus) stress, the transcriptome data were used. Heatmaps were visualized with TB tools. 65 4.6.3. qRT-PCR In order to explore the role of RgHDZ genes in responding to biotic and abiotic stresses, 11 differentially expressed genes were selected from different branches of each subfamily of R. glutinosa HD-Zip phylogenetic tree. The relative expression of RgHDZ genes in the tissue cultured seedlings of R. glutinosa under various abiotic and biotic stresses treatments were analyzed by qRT-PCR. Significant up-and down-regulated genes were determined as p < .05 using T test.
The total RNA of R. glutinosa was extracted with TRIzon (Beijing ComWin Biotech). BeyoRTTM III First Strand cDNA Synthesis Kit (Beyotime Biotechnology, Shanghai, China) was used to synthesize the first cDNA strand, qRT-PCR primers of RgHDZ and a reference gene TIP41 are listed in Supplementary  Table 3. The specific methods for qRT-PCR were performed according to Wang et al. 66

Conclusions
We identified 92 HD-Zip family genes from R. glutinosa. After bioinformatics and expression pattern analysis, we found that the HD-Zip family genes of R. glutinosa may not only regulate the growth and development of plants, but also participate in the biotic and abiotic stress process. The research showed that RgHDZ may respond to various stresses by down (or up)regulating its expression. Since, the entire growth and development stage of R. glutinosa, especially during the root expansion period, is vulnerable to various biotic and abiotic stresses, which greatly limits its yield and quality. This study analyzed the expression patterns of HD-Zip family genes in various coercions. Provide the theoretical basis for the cultivation of R. glutinosa resistant varieties.